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AN INVERSE -METHOD SOLUTION FOR RADIATING, NON ADIABATIC , 
EQUILIBRIUM INVISCID FLOW OVER A BLUNT BODY 

By Ralph A. Falanga and Edward M. Sullivan 
Langley Research Center 

SUMMARY 

An inverse flow-field program capable of making calculations with an equilibrium 
air gas model and a realistic radiation model has been developed. The program is capa- 
ble of calculating the inviscid flow field on a blunt body over a large range of conditions. 
The primary restrictions are that the assumed shock be analytic in shape and that the 
postshock temperature be between 10 000° K and 15 000° K. 

Comparisons are presented which show substantial agreement between the results 
from the present program, results from a direct-method program, and results from a 
time-asymptotic-method program for a nonadiabatic flow field. The inverse method can 
be operated with 8 to 24 rays and 5, 7, 9, or 11 radiative -flux calculations along each ray. 
A brief study has shown that a nominal set of 14 rays and 9 radiative flux calculations 
will produce accurate results over the entire range of applicability of the program. 

INTRODUCTION 

In the past several years a great deal of effort has gone into the development of 
computer programs which could be used to analyze the flow in the subsonic region behind 
the bow shock on a blunt body at hypersonic speeds. The early work in this area has been 
well documented, and the fundamentals can be found in a standard text such as that of 
Hayes and Probstein (ref. 1). In particular, three techniques have been developed and 
used extensively. The techniques, with an example of a working program for each, are: 
the direct method (Garrett, Suttles, and Perkins, ref. 2), the time-dependent method 
(Moretti and Abbett, ref. 3), and the inverse method (Lomax and Inouye, ref. 4). The 
solutions presented in these references are for inviscid adiabatic flows. (Ref. 2 carries 
a radiative flux term through the program development but never uses it in the 
calculations.) 

More recent work has been aimed at exploring the problems associated with entry 
velocities (10 km/sec or higher) at which the nonadiabatic effects become important. For 
stagnation-region calculations, a large number of specialized solutions have been devel- 
oped. In terms of the basic techniques mentioned previously, nonadiabatic effects have 



been included in the direct method by the work of Suttles (ref. 5) and in the time -dependent 
method by the work of Callis (ref. 6). The inverse method has not been as well developed. 
Cheng and Vincenti (ref. 7) used the classical inverse method to study radiation effects on 
various flow-field parameters. The analysis was, however, simplified by the assumptions 
of a perfect gas with a gray-gas radiation model. Olstad (ref. 8) used the simplified 
inverse method of Maslen (ref. 9) and an approximate nongray radiation model for equi- 
librium air. 

The present report presents the results of an effort to develop an inverse flow- 
field program which eliminates the approximations mentioned. The solution technique is 
the classical inverse technique described in reference 1. The gas model is equilibrium 
air, and the radiation model is the reasonably detailed model of reference 10, the same 
model as used in reference 5. This effort was undertaken primarily to provide solutions 
which could be used to test the validity of the more approximate programs, such as that 
of reference 5. The authors believe, however, that the effort was justified in its own right 
to fill the obvious gap in the inverse-program development. 

SYMBOLS 


Primed symbols (’) indicate dimensional quantities. Since it is convenient for solu- 
tion of the equations herein to have them in nondimensional form, the following charac- 
teristic quantities were used to nondimensionalize the individual terms: for lengths, 
the shock-wave radius Rg along the axis of symmetry; for velocity components, the free- 
stream velocity U'oo; for density, the free-stream density p'^; for pressure, twice the 
free-stream dynamic pressure p^U^; for enthalpy, the square of the free-stream veloc- 

O 

ity U’^; and for the radiation flux, the product of free-stream density and the cube of the 
free-stream velocity p^U^. 




H total enthalpy 


h static enthalpy 


J number of equally spaced points along each ray at which radiative flux is 

calculated 


K 


factor in governing equations, K = 1 - Qy 
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Q 


q R ,q R 
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T’ 

u,v 

V,U 


y 

z 
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V 
P 
<P 
w 

V 


static pressure 

curvature of reference surface (shock wave) 

radiative heat -flux vector and its magnitude, respectively 

body radius of curvature 

radius of curvature of bow shock at x = 0, Rg ~ ^R^ + 6^ 
radial distance from axis of symmetry of body (see fig. 1) 
temperature 

velocity components in x- and y-directions, respectively 
total velocity vector and its magnitude, respectively 
coordinate measured along reference surface (shock wave) 
a distance measured parallel to x-coordinate and downstream from it 
coordinate measured normal to reference surface (shock wave) 
axial coordinate from shock wave along axis of symmetry (see fig. 1) 
shock standoff distance 
transformed y-coordinate, rj = 1 - ^ 
density 

azimuth angle (see eq. (A7)) 
shock-wave angle (see fig. 1) 
del operator 



Subscripts: 


s condition at shock 

°° free -stream conditions 

o conditions along stagnation streamline 

SL sea-level conditions 


ANALYSIS 
Governing Equations 

The governing equations which describe an inviscid, nonconducting, radiating gas 


flow can be 

expressed in vector form as follows: 


Continuity 


p’(v . Vj + V • Vp' = 0 

(1) 

Momentum 


(y' . v)v = - -i- Vp’ 

(2) 

Energy 

pt 

/ f \ 



p’VV • VH ) = -v . 

(3) 

where 


TT t2 

H' = h’ + 

2 

(4) 


In the inverse method of solution, the shape of the bow shock is given, and the shape 
of the body and the details of the flow in the shock layer are unknown. A curvilinear coor- 
dinate system is used in the present analysis, and thus the generalization of equations (1) 
to (4) is removed by relating these equations to a reference surface, the shock wave as 
depicted in figure 1. The coordinate x is defined on the surface as distance along the 
surface from the system symmetry axis. The coordinate y is the distance from the 
surface along a normal to the shock, defined as positive in the direction toward the body 
from the shock. The curvature of the reference surface is denoted by Q(x) defined as 
positive as shown. 

The transposed equation (1) in terms of the shock-oriented coordinate system with 
dimensionless variables is written as follows: 
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Continuity 


a(pur) f a(pvrK) _ Q 

ax ay 

where K = 1 - Qy. 

Similarly, the momentum and energy equations can also be written as follows: 
x-momentum 


( 5 ) 


jfc +v J*i-Quv = --Li£ 
k ax ay k pK ax 


y-momentum 


av = .l^E 

K ax ay K pay 


Energy 
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(6a) 


(6b) 


(7) 


u 2 + v 2 

Note that replacing H with h + in equation (7) gives 


ah 2 9u dv Tr 3h it 8q „ 2 8v 
u — — + u^ — — + uv — — i- Kv — + Kuv — + Kv • 
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(8) 


The shock layer is treated as a one-dimensional gas slab for purposes of radiative 

9(q R r ) 

transport calculations. Thus, the term ■ in equation (8) is assumed to be negligi- 

ble. A detailed development of the form of these equations as used in the numerical solu- 
tion is given in the appendix. 

The gas model used is equilibrium air. Thermodynamic properties are obtained by 
using the curve-fit technique described in reference 2. 


Radiation Calculations 

In order to find a numerical solution to equations (5), (6), and (8), the radiative -flux 
9q R Kr) 

divergence term - — — - must be evaluated. Furthermore, one of the objectives of this 

program development is to use a reasonably detailed radiation model and avoid the use of 
gray gas or similar simplifications. The radiation model which was selected as being 
reasonably detailed but still compatible with the need for moderate computing times was 
developed for a radiation computer code called RATRAP (ref. 10). Local thermodynamic 
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and chemical equilibrium are assumed for the code, and a distribution of two thermody- 
namic variables must be known for the slab. (Pressure p' and enthalpy h' are used 
in the present analysis.) The radiation model for the RATRAP code is based on a non- 
uniform slab of high-temperature (T' ^ 10 000° K) air plasma. RATRAP, as noted, does 
not use a gray-gas approximation but takes into account the frequency dependence of the 
radiation. Included in the model are self-absorption in the gas, atomic-line radiation, 
and the important continuum radiation processes. The radiation processes included are 
the continuum, line emission, and absorption for atomic N and O, the ions N+ and 0+, 
electrons, and the O 2 Schumann-Runge continuum and N 2 Birge-Hopfield absorption band, 
This model was used without modification to calculate q^ wherever it appears in the 
present program. 


Numerical Procedures 

The numerical approach to this analysis is the inverse method, and the governing 
equations have been written in a form suitable for solution with this approach. The use 
of the inverse method for adiabatic flows has been well documented. (See, for example, 
refs. 1,4, and 11.) A similar numerical approach is adopted for the present program, 
which is for nonadiabatic flows. 

Briefly, the method consists of assuming a shock shape and a number of rays and 
their locations (fig. 1) and from the Rankine-Hugoniot relationships calculating the flow 
variables on each ray immediately behind the shock. The derivatives along the longitu- 
dinal coordinate (x-direction) can now be approximated by means of several mathematical 
schemes. For example, Van Dyke (ref. 12) used a series truncation method, whereas 
Marrone (ref. 11) used a seven-point polynomial fit to the calculated values. The latter 
approach is adopted for the present work (although the option exists of using a seven-point 
Lagrangian slope formula). Also, as an aid in delaying any tendency toward numerical 
instability, the flow variables are smoothed before taking derivatives in the x-direction. 
Smoothing was accomplished by using a least-square polynomial fit to each of the calcu- 
lated flow variables. This method had the beneficial effect of allowing larger integration 
step sizes with the result that overall computing time was reduced. Substitution of the 
evaluated longitudinal derivatives into the governing partial differential equations leads to 
a set of ordinary differential equations with the normal coordinate (y-direction) as the 
independent variable. These ordinary differential equations represent a set of equations 
which can be solved numerically (for example, by using the fourth-order Runge-Kutta 
integration technique) in the sequence given in the appendix. 

For an adiabatic flow-field solution by using the inverse method, the entire flow field 
is calculated for a given shock shape in a single pass. The presence of the radiative -flux 
term in the energy equation makes it necessary to resort to an iterative procedure in the 
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present program. The problem is begun with initial estimated distributions of pressure 
and static enthalpy along each ray, together with the length of the ray (i.e., the local 
standoff distance). This information is used as initial input to RATRAP, and initial 
values of the net radiative flux are calculated at J points along each ray. The local 
product Krq^ is formed at each point and derivatives in the y-direction are found by 
using a standard Lagrangian n-point differentiation formula. Since the tangent-slab 
approximation is being used, the resulting derivative is the flux divergence term required 
in equations (7) and (8). The solution proceeds until the "body point" is found for all 
except the last six rays. The resulting enthalpy profiles, pressure profiles, and standoff 
distances are then used as initial values, and the problem is repeated until a predeter- 
mined convergence test is satisfied, at which time the solution is considered to be com- 
plete. The convergence test is based on the static enthalpy at selected discrete points in 
the shock layer. Static enthalpy was selected for the convergence tests because it is one 
of the most sensitive parameters in the radiating gas flow calculations. The allowable 
error in h used for the enthalpy convergence test was ±0.006. Some calculations were 
made with the test value set at ±0.001, but the number of iterations greatly increased 
without any significant changes in the final results. No case has been found in which the 
iteration technique failed to converge to a solution. 

Since testing for convergence at every point in the flow field is not feasible, the test 
points were selected to be uniformly distributed throughout the flow field but included 
critical points where the enthalpy gradients are greatest. These test points do not, how- 
ever, include the last calculated y-station (body point) on any ray. It was originally 
planned to include these points, but the numerical problems associated with the mass bal- 
ance, which are discussed subsequently, precluded their use. 

An effort was made to establish the influence of the initial estimated profiles on the 
rate of convergence. The solution was found to converge for almost any reasonable esti- 
mate, but starting from the adiabatic profiles added one extra iteration. From experience 
it was found that the following procedure for establishing initial values gave consistently 
good results. For the pressure profile, the pressure is assumed to be invariant along 
each ray at a value obtained for the local oblique shock in a perfect gas. This value can 
be obtained from reference 13. For the standoff distance, the dimensionless standoff 
distance is assumed to vary with x according to the expression 

6 = 0.0392 + 0.022x (9) 

Finally, for the enthalpy profile, the enthalpy is assumed to vary with 6 so that 
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These equations, which were derived from adiabatic results and a linear approximation to 
the nonadiabatic results of reference 5, have been used extensively for both strongly and 
weakly radiating flows, and no case has been found in which convergence was not rapid. 

Three numerical procedure problems, which require some discussion, were encoun- 
tered. These problems are enthalpy profile damping, ray dropping, and satisfying the 
mass balance. 

Enthalpy profile d amping.- Initial calculations with the program showed that the 
solution did not converge in any reasonable number of iterations. The problem stemmed 
from the use of the updated enthalpy profile in each iteration without modification, and in 
some cases this use was sufficient to cause oscillation. This problem was handled by 
incorporating the following expression for updating enthalpies: 

h = h g - G(h g - h p ) 

where 

G damping coefficient 

hg enthalpy currently generated 

hp enthalpy previously generated 

Several values of G were tried and it was found that a value of G = 0.35 gave good 
results. This value has been incorporated into the program. The damping equation is 
not used until the start of the second iteration. Generally though, only two iterations are 
required after damping is introduced. 

Ray dropping. - As previously noted, the seven-point polynomial-fit technique (ref. 11) 
was used when taking derivatives in the x-direction. This technique was selected because 
of the excellent results obtained from its application to adiabatic flow, which does not 
require an iterative solution. The technique, as noted in reference 11, calculates deriva- 
tives until there are six rays which have not "reached the body" (i.e., satisfied the arbi- 
trary mass-conservation criterion). For an adiabatic flow, these last six rays are then 
discarded and the solution is considered to be complete. 

If carried over intact to the nonadiabatic flow, this numerical scheme would require 
that six rays be discarded at the end of each iteration. This procedure would seriously 
limit the number of iterations (maximum of three for the present program) and the 
x-distance covered by any solution. The calculated values from these rays are used only 
to provide updated p' , h' , and 6' solutions for RATRAP. To circumvent this difficulty, 
the initial total number of rays was retained, a procedure equivalent to assuming that the 
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last calculated values on each of these rays represented results at the body station. In 
the final solution the results from the last six rays are ignored since the calculations 
always show that they have not reached the body. The calculations do, however, show 
that the profiles along these rays are not so distorted that their use will influence the 
answers for the rays which do reach the body. 

Satisfying t he mass balance .- Since the inverse method starts with a given shock 
shape, one of the problems which is encountered in its use is establishing when the body 
is reached. The normal approach is to calculate the mass flow crossing the shock wave 
between the stagnation streamline and the ray in question and balance this against the 
mass flow crossing the ray as the flow is integrated in the y-direction. Reference 11 
shows that computing time becomes excessive as the mass-balance ratio becomes very 
close to 1.0, and the body shape associated with a ratio of 0.98 is indistinguishable from 
the body associated with a ratio of 0.995. Thus, the mass-balance test for a body point 
was initially set at 0.98 for the present program. It was soon discovered that because 
of the nonadiabatic nature of the problem, oscillations appeared in the enthalpy distribu- 
tion around the body. Furthermore, in some cases these oscillations tended to diverge 
with succeeding iterations so that it was not possible to obtain a solution. These oscilla- 
tions were not present in the enthalpy distributions for a mass-balance ratio of 0.96. 
Therefore, the mass-balance ratio requirement was lowered from 0.98 to 0.96 with the 
result that the oscillations did not appear and solutions were found for all cases. 

The influence of this change in mass-balance ratio requirement was investigated 
in the following manner. A series of solutions was calculated for several flight condi- 
tions with mass-balance ratio requirement values of 0.95, 0.96, 0.97, and where possible 
0.98. The stagnation-streamline (ray 1) results of these solutions were plotted as stand- 
off distance as a function of mass -balance ratio, and the resulting curves were linearly 
extrapolated to mass -balance ratios of 1.0. It was found that the standoff distance with a 
0.96 mass-balance ratio was 90.4 (±2.5) percent of the value for a mass-balance ratio 
of 1.0. This extrapolation technique has been applied to all values of 6 presented in 
this report except those for adiabatic flows for which the results with a 0.98 mass-balance 
ratio were considered adequate. Comparisons of enthalpy profiles from the present 
method with those from other methods, which will be shown subsequently, indicate that 
the present method is suitable for finding 6. 

In some cases, the solution did not oscillate but it did contain irregularities in the 
enthalpy distribution around the body at the last calculated y-station. These irregular- 
ities were traced to the seven-point numerical scheme used for taking derivatives in the 
x-direction. This scheme slides along the rays considering the first seven rays for which 
the body point has not been reached. When a ray reaches the body, it is removed from the 
calculation, and another ray is picked up to take its place. When the solution proceeds in 
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an orderly fashion so that rays drop out of the calculation one at a time, the numerical 
scheme produces a very regular enthalpy profile. However, in some cases the mass- 
balance ratio requirement is met on several rays during the same Runge-Kutta integra- 
tion step. When this occurs all these rays are removed at once, and the numerical 
scheme is forced to jump several rays; derivatives, which after one integration step are 
calculated by central difference formulas, are calculated by end-point formulas after the 
next integration step. It would be possible to overcome this difficulty by reducing the 
integration step size, but this has not been done because the error introduced was not 
sufficient to justify the additional computing time. For the static enthalpy profile, how- 
ever, the error may be sufficiently large to preclude the use of the last calculated 
y-station on each ray when calculating values for the convergence test. Consequently, 
these points have not been used when testing for convergence. 

Program Restrictions 

The program, as written, has several restrictions which must be noted. First, the 
assumed shock must be analytic in shape. Experience has shown that if the shock is not 
analytic, the program will not function properly. This restriction need not impair the 
usefulness of the program since many bodies can be found for which the shock approxi- 
mately meets this requirement. All the results given herein were calculated for the case 
of a catenary shock which produces a nearly spherical body. Second, the solution is con- 
fined to the subsonic part of the gas cap. Extension to the supersonic region would require 
conversion to a method-of-characteristics solution such as used in reference 4. 

Finally, as written, the program is restricted in the number of rays and radiative 
flux points J which can be used. The number of rays must be greater than or equal to 
8 and less than or equal to 24. The value of J must be 5, 7, 9, or 11. These limitations 
are not inherent in the inverse method but are the result of programing and machine stor- 
age requirements. 


RESULTS AND DISCUSSION 

The inverse flow-field program already described has been used with a catenary 
shock shape to calculate flow fields for a number of cases with a wide range of velocities, 
densities, and body sizes (table 1). In all cases the catenary shock shape produced a 
nearly spherical body. These cases, in general, cover flight conditions for which the 
postshock equilibrium temperature is between 10 000° K and 15 000° K and hence is within 
the range of applicability of both RATRAP and the air correlations of reference 2. The 
effects of varying the number of rays in the flow field and the number of points along each 
ray at which the radiation heat fluxes were obtained from RATRAP have also been studied 
to investigate the influence of the finite -difference mesh size on the calculated results and 
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the influence of the number of RATRAP calculation points on the accuracy of the flux 
divergence calculations. These cases are listed in table 2. In tables 1 and 2 the body 
radius was determined by subtracting the final calculated standoff distance 6* from 

the input stagnation-streamline shock radius Rg. Since the body radius is determined in 
this manner, it is not surprising that slight differences appear between the calculated 
R^ and the reference values in tables 1 and 2. 

Nonradiating Solutions 

In order to verify that the program was functioning properly, a nonradiating test 
case was conducted for comparison with the results of reference 2, in which the integral 
method was compared with the inverse method of Lomax and Inouye (ref. 4) . The com- 
parisons are shown in figure 2. As can be seen, the results from the three programs 
are generally in agreement. However, the results from the two inverse programs are in 
closer agreement, except for the density profiles close to the stagnation streamline. The 
differences that existed there, approximately 3 to 4 percent, are apparently the result of 
the different thermodynamic models used. The present program and reference 2 used the 
same air properties model whereas reference 4 used a different model. (Ref. 2 presents 
data which show that the air model of ref. 4 may be slightly more precise than the model 
used for the present work.) 


Radiating Solutions 

The influence of the nonadiabatic effect of radiation cooling, as well as the effects 
of body size and velocity, has been studied. The results are shown in figures 3 to 7. 
Figure 3 compares the nonadiabatic and adiabatic results calculated for a large body 
= 343.80 cm) at a high velocity = 14.55 km/sec). Various flow properties are 
plotted as a function of dimensionless shock-layer thickness for two rays, x = 0.009 and 
x = 0.254. The ray at x = 0.009 is the first ray in the mesh and will be referred to as 
ray 1. Since this ray is very close to the stagnation streamline, the calculated values for 
this ray have been assumed to be stagnation-streamline values. The ray at x = 0.254 is 
typical of a ray far from ray 1. Note that the pressure profile (fig. 3(c)) is not materially 
changed by the nonadiabatic effects, but all the other flow properties experience significant 
changes. As figure 3(e) shows, the standoff distance is decreased by approximately 
25 percent, but the resulting body remains nearly spherical. 

Figures 4 and 5 show the influence of body size on the ray 1 (stagnation streamline) 
enthalpy profiles for two different flight conditions. It can be seen that the nonadiabatic 
effects decrease with decreasing body size. The influence of the free-stream velocity is 
shown in figure 6, where the influence of body size and density have been eliminated. The 
decrease in nonadiabatic effects with decreasing velocity is evident. 
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Figure 7 shows a plot of enthalpy at the last calculated y-station for ray 1 as a 
function of free-stream velocity for two extreme values of free-stream density and for 
a large body. Note that the enthalpy extrapolates to the adiabatic value (0.5) for both 
densities at Ulc = 9 km/sec. From the results shown in figures 3 to 7, it is concluded 
that an adiabatic solution is always justified for flight velocities of 9 km/sec or less and 
for bodies with of about 3 cm or less. For flight velocities of 10 km/sec or higher 

and bodies with larger than about 3 cm, the nonadiabatic effects become significant, 

and an adiabatic solution could result in serious error. 

Radiating Solution Comparisons 

As previously indicated, one of the objectives of the present program development 
was to provide an independent check on other flow-field programs, particularly the method 
of integral relations or direct method (ref. 5). Comparisons have been made with the 
direct method (ref. 5) which uses the same thermodynamic and radiation models and with 
the time -dependent method from reference 6. Results from these comparisons are shown 
in figures 8 to 11. 

Figure 8 shows a comparison between the present results and those from refer- 
ence 5. The figure shows the flux distributions around the body for one adiabatic case 
R^ = 342.7 cm (ref. 5) and two nonadiabatic cases. The adiabatic case shows excellent 
agreement in the stagnation region (0 £ x 2 0.15), but shows a continually increasing diver- 
gence for x = 0.15. At x = 0.465 the difference is approximately 15 percent, based on 
the inverse solution. A comparison of the details of the two sets of results shows that in 
the stagnation region the standoff distance, pressure profile, and static enthalpy profile 
are identical, but at the larger values of x, differences appear in the standoff distance, 
shock curvature, pressure profile, and enthalpy profile. The flux distributions shown 
represent the cumulative effect of all these variations so that the spread shown is indica- 
tive of the comparative differences of the two methods for an adiabatic case. 

The comparative nonadiabatic flux distributions for two body sizes are also shown 
in figure 8. For 0 S x S 0.254, the results generally agree within 5 percent. For 
x > 0.254, the inverse -method results show the influence of modifying the ray-dropping 
procedure. The results from the last six rays (broken symbols in fig. 8) are known to be 
in error since the mass flow is never balanced (but never more than 10 percent out of 
balance) on these rays and the body point and flux profiles are never completely deter- 
mined. This kind of check shows that for nonadiabatic flow fields, either of the programs 
shown will produce radiative -flux values of comparable accuracy. 

Stagnation-streamline (ray 1) profiles of total enthalpy, density, temperature, and 
pressure for the nonadiabatic cases of figure 8 are shown in figure 9. The general level 
of agreement is good for both cases but definitely much better for the smaller body 
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(rJj = 34.27 cm). This effect is due to the previously noted fact that as body size 
decreases, the stagnation-streamline profiles approach the adiabatic values, which are 
in good agreement for the results from both programs. A comparison of the curves for 
the larger body = 342.7 cm) shows that as expected, the direct method does not give 
as much detail as the inverse method. Figure 9 also emphasizes that the inverse non- 
adiabatic calculations do not reach the body station. The last calculated point is at 
rj = 0.095 which is consistent with a mass -balance ratio of 0.96. 

In figure 10 are shown results taken from reference 6 for nondimensional enthalpy 
profiles along the stagnation streamline and those computed from the present program. 
Reference 6 uses a second-order time-asymptotic technique for computing a nongray, 
nonadiabatic absorbing flow field around a blunt body. The radiation model used in refer- 
ence 6 was the three -step model approximation. The results taken from reference 6 for 
comparison were based on a body nose radius of 2 meters and free -stream velocities of 
10 km/sec and 14 km/sec at a density ratio of 10"^ (equivalent altitude of 48.5 km). These 
results compared favorably with those computed from the present program for a velocity 
of 10 km/sec. (For this velocity, the radiation effects appear to be small but still non- 
adiabatic.) Also, for 14 km/sec where radiation effects are large, the two programs give 
results which differ at most by 10 percent. Some insight into this difference may be 
gained by examining the results for the low velocity. Here the radiation effects were 
small, and with the results in agreement, the conclusion can be drawn for nearly adiabatic 
flows that these methods for computing flow fields give comparable stagnation-streamline 
solutions. Hence, for the high velocity, the two remaining factors which can contribute to 
the 10-percent difference in results are the influence of strong radiation in the flow-field 
solution methods and the radiation models used. No attempt was made to isolate which of 
these two factors was dominant. 

Figure 11 shows results from three different flow-field computational techniques, 
the present inverse method, the direct method, and the time -asymptotic method. The 
effect of radiation on the enthalpy profile along the stagnation streamline is shown here 
for two body sizes, Rj, = 34.27 cm and R{j = 342.7 cm. The body sizes and flow con- 
ditions were the original ones used in reference 5. The direct method and inverse method 
used the same thermodynamic and radiation models while the time -asymptotic technique 
used the three-step radiation model and thermodynamic model described in reference 6. 
The results from the direct method and inverse method are in agreement for the small 
body (r{j = 34.27 cm) whereas for the large body (Rjj = 342.7 cm), the time -asymptotic 
scheme and inverse-method results agreed more favorably and displayed characteristi- 
cally similar enthalpy profiles. 

A final comparison with published results is shown in figure 12, where the flux and 
temperature profiles calculated by the present program are compared with the results 



from reference 14. Reference 14 considers a stagnation streamline only, integrates from 
the body to the shock, uses its own radiation and thermodynamic models, retains viscous 
effects in the momentum equation, and ignores the influence of curvature on the pressure 
gradient. A brief study showed that the dominant factor in the differences in figure 12 
is the continuum radiation model used in reference 14. This model uses a hydrogenic 
approximation which gives questionable absorption coefficients in certain spectral regions 
and distorts the flux profiles across the shock layer. The difference shown is similar 
to that shown in reference 15, which uses the inviscid flow-field calculations of refer- 
ence 5. In view of the agreement already shown (figs. 8,9, and 11) between the present 
method and reference 5, the differences shown in figure 12 are not surprising. 

Program Characteristics 

The results shown were obtained with the present program operating with what is 
considered a nominal mesh configuration, 14 rays spaced Ax = 0.0351 apart and 9 radia- 
tive flux calculations performed on each ray. The final part of the program development 
was devoted to an investigation of the rate of convergence and the influence of changing 
the number of rays, the ray spacing, and the number of radiative flux points on the final 
results. In all cases, the first ray was located at x = 0.009 from the axis of symmetry. 
The cases considered are given in table 2. 

The rate of convergence of the solution for a mesh with 14 rays and the radiative 
flux calculated at 9 points along each ray is shown in figure 13. Here the stagnation 
enthalpy profiles for the initial inputs through the fourth and final solution are shown. As 
illustrated by this plot, the solutions converge rapidly; convergence and the final solution 
are realized after the enthalpy accuracy criterion (±0.006) is satisfied. Note that for this 
case in the area of the body station the enthalpy convergence test is not satisfied between 
the third and fourth solutions because of the numerical procedural problem of forming the 
x-derivatives. The case shown took 8 minutes on the Control Data Series 6000 computer 
system. Computing times, however, are strongly dependent on the body size, free-stream 
velocity and density, the number of rays, and the number of points J along each ray at 
which the radiative heat flux is obtained from RATRAP. 

The effect of varying the number of rays from 10, Ax = 0.0507 to 20, Ax = 0.0240, 
for a fixed number of radiation heat -flux points along each ray is shown in figure 14 for 
a body radius of 34.44 cm. The effect on the stagnation enthalpy profile is negligible; any 
differences which may exist are within the order of accuracy for h (±0.006). These 
results are typical for even those profiles away from the stagnation streamline. Also, 
varying J from 5 to 11 along each ray and keeping the number of rays fixed at 14 for 
the same conditions above indicate negligible effects on the stagnation enthalpy profile , 
as shown in figure 15. However, the same conclusions cannot be drawn when the body 
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size is increased by approximately an order of magnitude. Here the effects are sizable 
in going from 10 rays and J = 5 to 14 rays and J = 9, as shown in figure 16. However, 
the differences become negligible when going from 14 rays and J = 9 to 20 rays and 
J = 11. On the basis of these results, it can be concluded that for general operation of 
the program, the nominal mesh configuration will produce good results over the entire 
range of applicability of the program. For smaller bodies or lower speeds, a much 
coarser mesh can be used to reduce computing time with no loss in accuracy. For the 
largest bodies and highest speeds considered, the gain in accuracy obtained by going to 
20 rays does not warrant the increased computer time. 

CONCLUDING REMARKS 

An inverse flow-field program capable of making calculations with an equilibrium 
air gas model and a realistic radiation model has been developed. The program is capa- 
ble of calculating the inviscid flow field on a blunt body over a very large range of con- 
ditions. The primary restrictions are that the assumed shock be analytic in shape and 
that the postshock temperature be between 10 000° K and 15 000° K. 

Results calculated by using the present inverse-method program have been com- 
pared with results obtained from a direct-method program and a time -asymptotic- method 
program. In general, the comparisons indicate substantial agreement among all three 
programs. However, the inverse method does retain more of the details of the flow than 
the direct method and therefore may be more valuable in certain instances. The degree 
of comparison between the inverse method, the time -asymptotic method of Callis (NASA 
TR R-299), and the stagnation-streamline solution of Rigdon, Dirling, and Thomas (NASA 
CR-1462) is clouded by variations in the thermodynamic and radiation models used. In 
general, the comparisons can be considered good. 

The inverse method can be operated with 8 to 24 rays and with 5, 7, 9, or 11 radia- 
tive flux calculations along each ray. A brief study has shown that a nominal set of 
14 rays and 9 radiative -flux calculations will produce accurate results over the entire 
range of applicability of the program. The use of more rays and/or more flux points 
results in no substantial increase in accuracy. The use of fewer rays and/or fewer flux 
points will save on computer time, but for larger bodies or higher flow velocities, a 
marked decrease in accuracy is noted. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., June 22, 1970. 
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APPENDIX 


DEVELOPMENT OF THE NUMERICALLY INTEGRATED 
GOVERNING EQUATIONS 

This appendix gives a detailed derivation of the governing differential equations 
(eqs. (5), (6), and (8)) and the procedure required to reduce them to a form suitable for 
sequential solution on a high-speed digital computer. The governing flow equations for 
an inviscid, nonconducting, radiating gas are given in vector form as 


Continuity 

p»(v . v ) + V • Vp' = 0 

(Al) 

Momentum 

( v • v) v = - -V v P ’ 

(A2) 


7 p 


Energy 

p’(v . VH’) = -V • q^ 

(A3) 

where 

11*2 

H* = h 

(A4) 


Primes in equations (Al) to (A4) refer to dimensional quantities. 

An elemental distance in the x-direction in terms of an elemental distance along the 
shock is (fig. 1) 

dx! = (l - Q’y’)dx’ (A5) 

and letting K’ = 1 - Q'y* reduces the elemental arc dxj to 

dx^ = K’ dx’ (A6) 

An element of arc ds' in three-dimensional coordinates, is 

ds'2 = K'^dx’^ + dy’^ + r'^dcp^ (A7) 

where cp is the azimuth angle. 

Thus, with this metric equation the relation between a general orthogonal coordinate 
system arid a specific shock-oriented one is established. Note that the coordinate system 
is curvilinear and orthogonal. 
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Then, the continuity equation (Al) in terms of the shock-oriented coordinate sys- 
tem with dimensionless variables and for an axisymmetrical condition can be written as 
follows: 


Continuity 

jfeSSi + *KP vr j9 = 0 (A8) 

ax ay 

Similarly, the momentum and energy equations can also be written as follows: 


x-momentum 


u 9u 9u Q 1 0p 

_ —Z + V — 1 - — uv = - - - 

K ax ay K pK ax 


y-momentum 


uav, v av,Q n 2 _ 1 9p 

K 8x 8y K pay 


Energy 


u 

9H 

T l = 

1 

[H 

V) 

K 

ax 

ay 

'pKr 

L ^ 


Note that replacing H with h + 


ay 


2 2 

- — in equation (All) gives 


(A9) 

(A10) 

(All) 


ax 


ax 


ax 


ay 


ay 


ay _ 

1 

[a ( 

w) 

9y 

"pr 

] dx 


9y J 


(A12) 


9/^0 r ) 

zSzJ. in equation (A12) 


ax 


The assumption of an infinite slab is used; thus, the term 
is assumed to be negligible. The energy equation is written 

u^l + u 2 |!i + uv-^ + Kv-ii + Kuv + Kv 2 -^ + = 0 (A13) 


ax ax ’ ' ax ay 

The equation of state is also utilized here 

P = P(h,p) 


9y 


ay pr dy 


(A14) 


These governing equations must now be manipulated so that the partial differential 
equations of the flow variables u, p, v, p, and h are obtained. Furthermore, the 
resulting equations should be written in a form suitable for sequential solution by using 
standard numerical techniques as mentioned previously. 

The term au/ay is obtained directly from the x-momentum equation (eq. (A9)). 


9u 1 (1 0p au 
ay Kv\p ax ax 



(A 15) 


Note that since the inverse method is being used, all derivatives with respect to x can 
be assumed known for the purpose of this derivation. 
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Next the y-momentum, energy, and state equations (eqs. (A10), (A13), and (A14)) 
are combined to obtain an expression for dv/dy in terms of dp/dy. The steps involved 
are as follows: 


From equation (A13), solve for 9v/9y 
1 


9 V 

ay 


Kv 2 L 


9 8a 8v 3h t/ 9h T/ 9u 1 

U Z TT- + UV — + Ut— + KV-7CT + Kuv — + — 


3 K Kr ) 

9y ' pr 9y 


ax ‘ ax ' “ ax 9y 

Differentiate equation (A14) with respect to y and solve for 9h/ ay 


(A16) 


8p _ j}p 8h + 8p 8p 
8y 8h 8y 8p 8y 


(A17) 


J_/9p _ 

ay “ c j \9y 2 ay j 


(A18) 


where Cj = 9p/9h and c 2 = 9p/9 p. The terms cj and c 2 are available from the 
correlation of equilibrium air properties to 15 000° K described in reference 10. Sub- 
stituting the expression for 9h/9y (eq. (A18)) into equation (A16) gives 


9v 

ay 


l 

Kv 2 


2 9u 9v 9h Kv /9p 9p\ , v „ 9u 

9x 9x ax Cj y9y 2 dyj 9y 


1 8 (' 

*R Kr ) 

pr 

ay 


(A19) 


Replacing 9p/9y in equation (A19) by its equivalence from the y-momentum equa- 
tion (A10) and collecting all terms in 9v/ 9y gives 


9v L _ p\ _ 1 [.2 au 9v ah _ puv dv _ Qpu 2 v _ Kvc 2 

ay\ C 1 J Kv 2 U ax + 9x + 9x c^ 9x Cj Cj 9y 


9u , l 8 ( <1 : 

+ KuV 3y + pr 9y 


R Kr ) 


(A20) 


Equation (A20) is the desired equation for 9v/9y in terms of dp/dy. Equations (A20) 
and (A8) now form a set of two equations in two unknowns which can be solved as follows: 

From equation (A8) 
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The expression in equation (A20) for 9v/8y is substituted into equation (A21) to give 


ap 

ay 



/u 3p 9u 
\p ax + ax 


. u 9r 
r dx 


+ 


t^-Q v 

r 9y 



X 



QPU^V + Kuy 9u + 

c x ay 


1 a (qR Kr ) 

pr 9y 


+ u 


2 9u 
9x 


(A22) 


Thus, once dp/ dy has been evaluated with equation (A22) , then 9v/9y can be evaluated 
with equation (A20) . Now all information is available to evaluate dp/ 9y from the 
y-momentum equation (eq. (A10)). Finally, ah/8y is evaluated from equation (A18). 


In summary, the final set of partial differential equations (eqs. (A15), (A22), (A20), 
(A10), and (A18)) used in the program, listed in the order in which they are evaluated, is 


9u _ 1 fl 9p 

dy Kv \P dx. 


+ u-g-Quv) 


9P 

9y 



+ Kv|£.q v 

r dy 


+ 



X 


uv 



dv . „2 ah _ Qpu 2 v 

ax ax ax c x 


+ Kuv 

9y 


1 9 ( q R Kr ) 

pr 9y 


dv 


Kv 2 / 1 - £ 


T , 8u 1 
+ Kuv — + 


9 8u 9v ah 
u“ + uv — + u — — 
ax ax ax 


3i 


( ( lR Kr ) 


ay pr ay 


iP = -pfH. il + v^ + Su^ 

ay P \K ax ay k ) 


ah = i_/a£. c 
ay " c i \ay 2 ay 


puv dy Qpu 2 v 
C 1 ax " c 1 


Kvc 2 a P 
c i ay 
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TABLE 1.- RANGE OF VELOCITIES, DENSITIES, AND BODY SIZES 
STUDIED INCLUDING COMPARATIVE CONDITIONS 


P °°/ PSL 

Uoo, 

km/sec 

cm 

Rjj, cm 

Remarks 

Inverse 

method 

Reference 

2.2064 X 10" 4 

15.25 

3.105 

31.9 

314.85 

2.971 

30.7 

304.4 

3.05 

30.5 

305 

Velocity and body size (by two 
orders of magnitude) varied at 
altitude of 61 km 

2.2064 X 10-4 

12.81 

314.85 

302.9 

305 

2.2064 X lO- 4 

10.37 

314.85 

3.105 

300.8 

2.961 

305 

3.05 

1.392 X lO' 2 

11.59 

3.105 

31.9 

314.85 

2.945 

30.30 

300.3 

3.05 

30.5 

305 

Velocity and body size (by two 
orders of magnitude) varied at 
altitude of 30.5 km 

1.392 X lO -2 

9.455 

31.9 

314.85 

30.4 

298.6 

30.5 

305 

1.84 X lO -4 

14.55 

356.16 

35.84 

343.8 

34.44 

342.7 

34.27 

Conditions of reference 5 

1.0 X lO' 3 

14 

10 

208 

209.5 

199.50 

199.50 

200 

200 

Conditions of reference 6 



TABLE 2.- IN FORMATION USED TO ESTABLISH FLOW- FIELD MESH SIZES 



T 

t 



Number of 


P °°/ P SL 

Voo, 

km/sec 

K b> 

cm 

Number 
of rays 

Spacing of 
rays, Ax 

radiative flux 
points along 
each ray, J 

Remarks 

1.84 X 10' 4 

14.55 

34.44 

20 

2.401 x 10" 2 

5 

Number of rays varied, 




10 

5.07 X lO -2 


total x-distance fixed 




14 

3.509 x lO -2 


at 0.464, and J = 5 

1.84 x 10' 4 

14.55 

34.44 

14 

3.509 x 10-2 

9 

Only J varied 


11 

1.84 x lO' 4 

14.55 

34.44 

10 

5.07 X lO" 2 

11 

Effect of minimum num- 







ber of rays used and 
maximum J 

1.84 x lO" 4 

14.55 

34.44 

20 

3.509 x lO -2 

5 

Increased x to 0.675 







for minimum J 

2.2064 X 10" 4 

15.25 

304.4 

10 

5.07 X lO -2 

5 





14 

3.509 x lO" 2 

1 

9 

Effect of changing num- 




20 

2.401 x lO' 2 

11 

ber of rays and J on 
large body and high 







2.2064 x 10" 4 

12.81 

302.9 

10 

5.07 X lO" 2 

5 

velocity and altitude 




14 

3.509 x lO" 2 

9 
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Figure 1.- Schematic of coordinate system and flow field. 








(a) Enthalpy. 


(b) Density. 


Figure 3 


x=.009 


x= . 254 


1.0 r 


p 



Adiabatic 

Non adiabatic 


J 1 i I 

.02 y .03 .04 .05 

(c) Pressure. 


Flow parameter profiles for adiabatic and nonadiabatic conditions. 

p /p_ T = 1.8k x 10-^j R-k = 3^3.80 cm. 


» = 1^.55 W‘ 




0 .01 .02 .03 .04 .05 0 .1 .2 .3 

y Z 

(d) Temperature. (e) Generated bodies for catenary 

shock shape. 


Figure 3.- Concluded. 
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Figure 5.- Variation of stagnation-streamline 
enthalpy profile with R^. U^, = 15.25 km/s 
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Figure 6.- Variation of stagnation-streamline enthalpy 
profile with = 305 (reference value); 

P c»/ P SL = 2>20 ^ X 10-^; altitude = 6l.O km. 
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Figure 7-- Effect of free-stream density and velocity on stagnation-point enthalpy 

R-k = 305 cm ( reference value ) . 
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Figure 8.- Comparison of radiative heat flux along body surfaces from present 
inverse program with results from the direct method. Hi = 1^-55 km/sec; 

P„/P S L = X - 84 X 10 " 4 • 
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Comparison of stagnation-streamline flow parameter profiles from the present inverse program with results 
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Figure 10.- Comparison of stagnation-streamline enthalpy 
profiles from inverse method with those from time- 
asymptotic method (reference 6). Poo/PgXj = 

R-/ = 200 cm (reference value). 


Figure 11.- Comparison of stagnation-streamline enthalpy 
profiles from present program with those from refer- 
ence 5 and reference 6. Uj, - 1^.55 km/sec; 

P m /P S L - ^ X 10 " • 
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